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ABSTRACT 


We have a system of N sendox machines (service centers) 
which can transmit or receive (but not simultaneously) photo- 
copies to and from each other. The stochastic process 
{X,,(t),t > O}, where X(t) is the number of messages under- 
going transmission in the system at epoch t, is selected to 
be the measure of system congestion. Under the assumptions 
of exponential inter-arrival and transmission times and that 
ie not transmitted immediately are lost, X(t) Ag a 
birth and death process. We call the model under these as- 
sumptions the "loss" model. It is shown that Xt) can be 
approximated by a diffusion process. The steady state re- 
sults of the diffusion model can be calculated with ease and 
are in close agreement with those of the birth and death 
model and the simulation "loss" model. A more realistic "no 
loss'' model is also developed where messages that are not 
transmitted immediately are allowed to queue up. It is shown, 
however, that this “nO loss" model can be quite accurately 


approximated by an appropriate "loss" model. 
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NOTATIONS 


N total number of machines (Service centers) in the sys- 
tem. (ASSumed to be an even number. ) 


M=N/2 maximum number of messages that can be transmitted by 


the system. 

arrival rate of messages at each machine. 

ort effective arrival rate of messages at each machine in 
the "loss model. 

Le departure rate of a message undergoing transmission. 


X(t) number of messages undergoing transmission at epoch t. 
S(t) Stochastic element associated with Xt). 

y(t) deterministic component of X(t). 

,(t) characteristic function of X(t). 

Py, Ct) characteristic function of Sy(t)- 

< Steady state mean of X(t). 

Steady state variance of X(t). 


Pos Steady state probabilities of X(t). 
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I. INTRODUCTION 


The main purpose of this thesis is to construct and com- 
pare various models which can be used to analyze the degree 
of congestion in a system of N sendox facsimile TROMASS 
sendox machines can both transmit and receive (but not simul- 
taneously) photocopies uSing standard telephone lines. The 
original problem arose out of a proposal to eran ra number 
of these machines in central London Mintech buildings for 
transmission of thin "immediate" messages between buildings. 
This problem was first studied by Coleman [Ref. 1] who in- 
vestigated the expected system congestion, expected cost per 
message and a possible priority system. 

One of the models in his study considered each machine 
individually as a single server system. It was assumed that 
only one machine was located in each building. Arrivals at 
each of the N machines occur from two sources: 1) requests 
to transmit messages to other machines and 2) incoming mes- 
sages from the other N-1 machines. These two arrival streams 
were merged into one and the standard queuing formula for a 
Single server queue was applied to obtain the expected wait- 
ing time of messages at an individual machine. However, a 
subsequent simulation study showed wide discrepancy between 
the theoretical and simulated expected waiting times. This 
occurred because the Single server model did not take into 
account the "interference' between machines. The model as- 


sumed that when a message is transmitted from one machine, 
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its destination (which could be any one of the other N-1 ma- 
chines) is always available to receive the message. In 
reality, its destination may be engaged in transmitting or 
receiving another message. The effect of such "interference" 
tends to increase the waiting time as compared to the result 
obtained from the single server queuing model. Any reason- 
able theoretical model, therefore, has to consider all N ma- 
chines simultaneously. It is very difficult, however, to 
build an.analytical model that can give a reasonable esti-~ 
mate of waiting time because it has to consider N queues all 


at once. 


neal 


II. DESCRIPTION OF THE "LOSS" MODEL 


The problem becomes much more tractable when we only con- 
sider modelling the system size i.e., the number of messages 
undergoing transmission. This number varies from 0 to M 
where M=[N/2], the largest integer less than or equal to N/2. 
Let Xt) be defined as the number of messages undergoing 
transmission at epoch t. The stochastic process {X(t ), 

t > 0} will be our measure of system congestion. 

A number of assumptions are made to facilitate construc- 
tion of theoretical models for the process X(t). We assume 
a Poisson arrival of messages at each machine. The arrival 
rate of messages requesting transmission at each machine is 
the same and is denoted by Me Moreover, we assume that each 
message undergoing transmission is finishing at an exponen- 
tial nee of uw. Furthermore we define p to be = and hence 
the departure rate (yu) of a message is also equal to 2pd. 
Finally, we assume that if a message fails to be transmitted 
immediately upon arrival (due to either the receiving or 
sending machine being busy), then the message is lost to the 
system. For this reason, we call the model based on the 
above assumptions the "loss" model. 

Two theoretical and two simulation models for the proc- 
ess X(t) will be described and discussed in the following 
sections. We are mainly interested in steady state behavior 
of the process X(t). In the first theoretical model, X(t) 


1s shown to be a birth and death process (thiS was done in 


a2 





{[Ref. 1]) whose steady state probabilities can be derived 
and. calculated numerically. Since Iglehart [{Ref. 3] had 
proved the convergence of certain birth and death processes 
to diffusion processes, we were led to consider a diffusion 
approximation for the process X(t), when M is large. This 
was the approach taken in the second theoretical model where 
it was shown that X(t ) can be approximated by a non-station- 
ary Ornstein-Uhlenbeck process. Finally, a simulation model 
was built to validate the theoretical models. The steady 
state results of the theoretical and the simulation models 
compared very favorably with one another. 

A second simulation model was developed for the case in 
which messages not transmitted immediately, queue up until 
they are finally transmitted. We call this model the ‘no 
loss" model. We are unable to construct a theoretical model 
for this case. However we will show that this "no loss" 
model can be quite accurately approximated by are diffusion 


"loss" model. 
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III. BIRTH AND DEATH "LOSS" MODEL 


In an unpublished paper [Ref. 1], Dr. Rodney Coleman 
modelled the number of messages undergoing transmission, 
X(t), as a birth and death process. For the sake of sim- 
plicity, we assume that N, the total number of machines in 
the system is even and hence M, the maximum number of mes- 
sages that can be carried by the system is given by N/2. 
When the process X(t) is in state jy, the conditional proba- 
Dilaty that X(t) will visit state jt+l in the next interval 
of length At is ee + o(At) and the conditional probability 
that X(t) will visit state j-1 is Bee + o(At), where a, 
and Ps are the infinitesimal birth and death rates respec- 
tively [Ref. 6]. MesSages can only enter the system at the 
(2M-2j3) free machines at a rate of (2M-2j3)A. When a message 


arrives at a free machine, its probability of finding an 


available receiving machine is a, therefore, 
— (2M-23 )(2M-2j-1)A 
Me (2M-1) (Sets) 


Since there are j messages, each finishing transmission at 


the rate of u, and thus 


We SW a asides (3.2) 
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Let es denote the steady state probability of X(t) being in 


state j3. Then Coleman [Ref. 1] showed that 


k 
ie... pau i: 
= K1C2M-2k)t = 
p, = EN Sh k=0,1,...M 
yk: [ge ee ee 
where Y = (oM-a) Vb * BpCeii-1) 


However, Since the caluclation for the PLS are carried 


out on a digital computer, a more efficient method computa- 


tionally is to use the following recursive relationships 


fRef. 6]. 
Peete ore Die k=1,2,...M C52) 
M , 
KZ0°K ~ 7 (2082) 


The computer program which calculates the steady state prob- 
abilities and also the mean and variance of X0t) at steady 


state is listed in Appendix D. 
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IV. DIFFUSION "LOSS" MODEL 


We follow Gaver's [Ref. 4] approach by approximating the 
process X(t) with a deterministic component y(t) and a 
stochastic element or noise S(t). Thus X(t) can be written 
as: 


X(t) = My(t) + JM S(t) (4.1) 


We show that as we let M go to infinity, S(t) converges 
to a non-stationary Ornstein-Uhlenbeck process, S(t) [Ref. 
2). Hence for large M, we can use the following formula to 


approximate X(t): 
X(t) = My(t) + YM S(t) (42) 


Because of our assumption of Poisson arrivals, a message 
will arrive at a machine with probability Adt + o(dt) in the 
time interval (t,ttdt). During this time interval, the 
probability of two or more messages arriving at the same ma- 
chine is o(dt). At epoch t the number of messages under- 
going transmission is X(t) and hence the number of unoccupied 
machines iS given by N-2X,,(t ) or 2M-2X,,(t). In the time in- 
terval (t,ttdt), effective arrival of messages can occur 
only at any of the 2M-2X),(t) unoccupied machines. Given that 
a message has arrived at an available machine, its probabil- 
ity of finding the receiving machine not in use is now 
(2M-2X,(t)-1)/(2M-1). Therefore, 


Prob{X ttdt )=X,,(t)+1|given X(t )} 


( 
M 
(2M-2X,,(t))(2M-2X, (t)-1)Adt (ws) 


J . ———————————————————— + 
OM—1 o(dt) 


16 


SS => == 





Similarly, under the assumption of exponential transmis- 
Sion times, each message undergoing transmission is finish- 
ing at the rate of uw or 2pA. At epoch t, messages will be 
finishing at a rate of 2X,,(t) pr. Therefore during the time 
interval (t,ttdt) the probability of a message completing 
its transmission is 2X,,(t ) prdt + o(dt) and the probability 
of two or more messages completing their transmission is 


o(dt). Therefore, 


Prob {X,,(ttdt) = X(t )-1 | given X(t) } 


(4.4) 
= 2X,,prdt tO at) 
Finally, 
Prob {X,,(ttdt) = X(t) | given X(t) 
(2M-2X (t)) (2M-2X,(t)-1) (a8), 
oe eo ee et pAdt + o( at) 


In Appendix A, equations (4.1) through (4.5) are used to 
derive the partial differential equation 
9 Z 
SY = _an{atp-2y(t)}o2% - 2aty2(t)-y(t)(2-p)+1}2y 
ot 36 2 
(4.6) 
Here y~(6,t) is the characteristic function of S(t), the lim- 
it of the noise process S46") as M> ©, The deterministic 


process y(t) must satisfy the ordinary differential equation 
y'(t) - 20{fy*(t) - (2t+p)y(t) - 1} = 0 (4.7) 


in order that the partial differential equation (4.6) be 
valid. The partial differential equation (4.6) is recognized 


as the transformed version of the forward differential 


Inf) 





equation for a non-stationary Ornstein-Uhlenbeck process. 
It is Known [Ref. 2] that the probability density function 
of the Ornstein-Uhlenbeck process S(t) must satisfy the 


Kolmogorov Forward Differential Equation 


2 
SP = 4 2) {a(s,t)p} ~  {8(s,t)p} (4.8) 


where £(s,t) and a(s,t) denote the infinitesimal mean and 


variance respectively, conditioned on S(t)=s. 


BCs. t) -2A{2+p-2y(t)}s (4.9) 


a(s,t) = 2r{y*(t)-(2-p)y(t)+1} (4.10) 


The first step toward determining the mean and variance 
of the process S(t) involves solving the first order differ- 
ential equation (4.7) to obtain y(t) explicitly, assuming 
the initial condition that y(Q)=0. The solution is derived 
in Appendix A and is 


—-4,wt 


1 + ke } (4.11) 


@ Awl 


y(t) = v-w { 
ek 


where v=1 + #p 


w? = (1 + 49)? - 1 


Ww-V 
wtyv 


Since the Ornstein-Uhlenbeck process is Gaussian, we may 
differentiate the Gaussian characteristic function 
exp [ida(t)-#067b(t)] with respect to 6 and t, substitute jy, 
ot oY back into equation (4.6) and equate the coefficients 


of i160 and 6? in order to obtain two first order differential 


equations for a(t) and b(t) which denote the mean and 
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variance of the process S(t) respectively. Assuming that 


S(0)=0, we can show that 
E[S(t)] = 0 t > O (4.12) 


and since 


X(t) = My(t) + YM S(t) 
we obtain 
E[X,(t)] = My(t) +t > 0 (4.13) 


where y(t) is given by equation (4.11) 

The solution of the first order differential equation for 
the Var {S(t)} involves integrals that are quite intractable 
and we are only able to obtain an approximation for large t. 


2 
Yo ~ ¥,€2-p)t1 


2 1c + f AN 
Var {oct 2(D*p-2y_) t large ec 
= lim = 1 1 2 
where Ving eas y(t) = 1+ 4p - V(1+4p)7-1 C4155 
Since Var {S(t)} = Var {X(t}, we obtain the following 


expressions for steady state mean and variance of X(t), de- 


noted by x, and oe respectively. 


x = E[X,,(t)] = My t large 
: m e (46) 
ae 
me Var {X(t )} 
Dey Cee al 
= M ] t large C4417) 


'Z(B¥ p= By) 


Detailed derivation of equations (4.12) and (4.14) is given 
in Appendix C. When t is large, X(t) is normally distri- 


buted with mean and variance given by equations (4.16) and 


Lg 





(4.17)respectively. Plots of x_/M and 7 /M against P are 
given in Figures 1 and 2 respectively. It is seen in Appen- 
dix C that y(t) reaches its steady state value of vies consid- 
erably: faster than Var {S(t)}, thus the steady state mean of 
XC t) is attained much faster than the steady state variance 
of X(t). 

The steady state probabilities for Xt) can be easily 
obtained from the standard normal table after x_ and o* are 
calculated using equations (4.16) and @.17). However, since 
we are uSing a continuous (normal) probability distribution 
to approximate a discrete distribution, we have to make the 


following "continuity corrections" 











SRS) 
Pos SS ab (4.18) 
2710 
‘ Ca). 
a rT cae 
Pp. = sla e S dx 1 alt ae |) col (4.19) 
me /2c6 
(x-x_)° 
(ee) os 
Pp, = ip = eo <6 dx (4.20) 
M-3 Y270. 


Another way to derive the parameters for the Ornstein- 
Uhlenbeck process S(t) is by defining the stochastic differ- 


ential equation for the process X(t). 


dX, ,(t) = By Cs, t)dt = Lt.) van (x, t)dt (4.21) 


where Z(t) is a purely random Gaussian process with zero 
mean and unit variance and B(x, t)dt and a(x, t)dt are the 


mean and variance of the increment dx, ¢t) in a small time 
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interval (t,ttdt). Recall that we have derived the follow- 
ing probabilities earlier 
Prob (dX, (t) = 1] given X(t) 3 


: {2M-2X,,(t)}{2M-2x, (t)-1} 


SM Aat + O(dt ) 


Prob {dX,,(t) = -1| given X(t )3 


= 2X,(t)pddt + o(dt) 


Prob {dX,,(t) = 0O| given X(t) 3 


{ 2M-2X,,(t )}{2M-2Xx,,(t)-1} 
= 1 = CAG O= 2%, (t) pdt + (dt) 


2M-1 
Therefore it is reasonable to define B,, (x, t)dt and a(x, t)dt 
as follows: 
By (x, t)dt a E[ dx, (t )] 


(4322) 
{2M-2X,,(t) }12M-2x, (t)-1} 
eG Ol 


QM-1 
a(x, t)dt = Var [dX,,(t ) ] 
(4523) 
{2M-2x,,(t)}{2M-2x,,(t)-1} 
aan em Adt+2X,(t )prdt 


Substituting 54, (t) as given by equation (4.1) in place of 
X(t), we derive a stochastic differential equation for 
Sy Ct) in Appendix B. Then by letting M > » and requiring 
the first order differential equation (4.7) for the deter- 
ministic process y(t) to hold, we obtain the stochastic 


differential equation 
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ds(t) = 2A{2y(t)-2-p}S(t)dt 


+ Z(t) V2A{y7(t)-y(t)(2-p)t+1}dat (4.24) 


This equation is the stochastic differential equation for an 
Ornstein-Uhlenbeck process which has the same infinitesimal 


mean and variance as those given in equations (4.9) and (4.10). 
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V. SIMULATION "LOSS" MODEL 


A. SIMULATION PROGRAM I 

Two methods of simulations are presented here for the 
"loss'' model. Simulation Program I makes use of the obvious 
approach of duplicating the operation of ae actual system. 
Messages arrive at each of the N machines at a Poisson rate 
A. These messages are either lost or undergo transmissions 
which are being completed at an individual rate of u. A 
brief description of the program logic is given in the next 
paragraph. 

Simulation Program I is written in FORTRAN. It operates 
simply by searching for the next event to occur. [In this 
model an event is either the arrival of a message at a ma- 
chine or the completion of a message undergoing transmission. 
Each event iS assigned a clock time which denotes the time 
of occurrence of that event. The next event is that whose 
clock time is the minimum of the clock times of all the fu- 
ture events. Having found this event, the program performs 
a series of tests to determine whether the conditions in the 
System are such that the potential event can take place. If 
so, the program moves to execute this event. For example 
the event may be the arrival of a message at machine i, then 
the necessary test is to check whether the sending machine 
1 and a randomly selected receiving machine j, j#i are both 
available. If so, the number of messages undergoing trans- 


mission increases by one. The basic loop is repeated as 
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many times as necessary to complete the entire Simulation. 
The flow chart and listing of this program are given in Ap- 


pendix D. 


B. SIMULATION PROGRAM II 
Another way to develop a Simulation program for the "loss" 

model is to make use of the fact that X(t) is a birth and 

death process. We know the process X(t) Stays in state i 

for an amount of time that is exponential with rate Aj TUS 

and given that it makes a transition from state i, it goes 

to either state i-1 or itl with probabilities are 


. Age 
z respectively where ds and Hs are as defined in section 


and 





AqtHa 
ie. 


Instead of monitoring the flow of messages through the 
system of N machines, Program II works by keeping track of 
the length of time that the process X(t) spends in each 
state and when X(t ) changes state. This program is also 
written in FORTRAN. It starts off with XK, C0) in a given 
initial state i. An exponential variate with mean ~~ 
: Agtiy 
is generated and this represents the length of time that 


X(t) spends in state i before jumping to state i-1l or itl. 


In order to decide whether X(t) next visits i-1 or itl, a 





HittOoumemandem vanrate p, O <p < 1, is generated. Ji p is 
less than x = Dene wt we viscite state afl next and if p 
AFH A M 
exceeds <= X,,(t) goes into state i-1. Once X,,(t) switches 
itu OM | 


to a new state j, we repeat the process of determining the 


length of X(t )'s Stay in state Jj) and the next state to be 
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visited. If X(t) is in state 0, it will next visit state 
1 with probability one. Similarly if X(t) is in state M, 
it will next visit state M-1 with probability one. The out- 
put of this program will be discussed together with that of 
Program I in the next section. The flow chart and listing 


of Program II are included in Appendix D. 


C. DISCUSSION OF OUTPUT OF PROGRAMS I AND II 

Both Simulation programs have five input parameters, 
namely the total number of machines (N), the duration of each 
Simulation run (TSTOP), the arrival rate of messages at each 
machine (A), the departure rate of transmitted messages (1) 
and the initial state i. Since we are mainly interested in 
steady state results, the primary output of both simulation 
programs is the probability distribution of the process 
X(t) when in steady state. 

One way to obtain these steady state probabilities is to 
let the simulation program run for t' minutes until we are 
reasonably sure that steady state has been achieved, and 
then record that value of X(t! )- Since each Simulation run 
produces only one realization of X(t"), this method of es- 
Eilnatving the wg proves to be very inefficient because of 
the large number of simulation runs required to produce a 
consistent estimate. 

Papenlctie Metnodn tO sccbamate the ee is to observe the 
time average probabilities of X(t ) when it is already in 
steady state. We choose a starting state 1 which is very 


close to the steady state mean MY, so that X(t) will reach 
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steady state shortly. The simulation model will have to run 
for a considerable amount of time in order to gather enough 
sample values. However, this results in a program that re- 
quires too much core storage. To overcome this problem, we 
made 200 simulation runs of TSTOP ninieee instead of a single 
run of 200 X TSTOP minutes. 

During each execution of the simulation program, 200 
Simulation runs of duration TSTOP minutes each are made. 
Let f. denote the frequency that the process X(t) visits 
state i. At the beginning of the program, all the f.'s are 
set equal to zero. During Eel Simulation run, the number 
of messages undergoing transmission is observed at discrete 
one minute intervals. If this number is k, then f ls incre- 
mented by one and all ee fa: j#k, remain unchanged. 
Thus, during the entire simulation program, X(t) is observed 
at 200 x TSTOP points, and the estimated steady state proba- 
bilities are given by: 


A a 


Pr ~ 900xTSTOP 


Initially, both programs are tested with the same input 
parameters to ensure that they are yielding results that are 
consistent with each other. From then on, only Program II 
is used to generate the various output that are used to val- 
idate the theoretical models, because time-wise, it is much 
more efficient than Program I. For example, with N=50, 
A=2.5 and w=10.0, the total execution time of Program I is 


about 18 times that of Program II. 
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Vi. SIMULATION "NO LOSS" MODEL 


In the previous sections, we have developed analytical 
and ssinculais en models for a system of N machines in which 
messages that are not transmitted immediately are lost. A 
more realistic assumption is to allow these messages that 
are not transmitted immediately to join a queue at the re- 
spective machines. We call a model based on this less strin- 
gent assumption a "no loss" model. This "no loss" model is 
not easily amenable to theoretical analysis, but it can be 
eaSily simulated by a program written in GPSS. The assump- 
tions of exponential distribution for the inter-arrival and 
inter-departure times of ‘messages could now be relaxed be- 
cause a Simulation program can easily handle any other prob- 
ability distribution. However, the previous assumption of 
exponential inter-arrival and transmission times is retained 
for comparison purposes. 

In the GPSS program, the Sendox machines are treated as 
facilities and the messages that arrive at the various ma- 
chines are treated as transactions. The program consists 
mainly of N closely similar segments. Each segment simulates 
Mac aclivicy Of a facility. The program instructions for 
each segment are almost identical except for the fact that 
Eaellevaciiatyehas a sdarierent facility number i, 1=1,2,...N. 

Arrivals of messages are simulated through the generation 
of random arrivals of transactions at an exponential rate Xi 


memcaciolneciwity. “A transaction created at facility i is 


NS, 





randomly assigned to a destination facility j, j#i, inde- 
pendently of other meen tone: Then facilities i and j 

are tested to see whether they are available simultaneously. 
If so, the transaction proceeds to seize both facilities for 
a length of time that is exponentially distributed with rate 
Hw. If not, the transaction will have to wait until both 
facilities i and j are simultaneously free. If another trans- 
few Tonle seecreatca while a previous transaction 1S either 
waiting for or undergoing service, the new transaction must 
queue up. Therefore we assume a first come first serve 

queue discipline. This procedure is followed at all N facil- 
ities. 

GPSS is preferred to FORTRAN in this model because the 
logic and structure of : GPSS program are inore Simple. Moyre- 
Aver GHEE statistics and steady state probabilities of 
X(t) can be gathered with little extra effort. The only 
Shortcoming is that this program cannot adapt easily to 
changes in the number of machines. Thus if we want to double 
the number of machines the card deck will almost double its 
size. This GPSS program was only written for N equal to 20, 
whereas N can go as high as 200 in Programs I and Il. A 


listing of this program is included in Appendix D. 
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VII. COMPARISON OF RESULTS 


According to the diffusion "loss" model, when X(t) is 
in steady state X(t) is approximately normally distributed 
with mean x, and variance oe given by equations (4.16) and 
(4.17). Moreover, using x_ and oi and the continuity correc- 
tion, given by equations (4.18), (4.19) and (4.20), the 
steady state probabilities can be obtained from the standard 
normal table. 

We are mainly interested in determining the accuracy of 
the diffusion "loss" model for various values of N and p in 
comparison with the simulation models. For ease of refer- 


ence equations (4.16) and (4.17) are repeated here. 


= aa 
X My. 5 V, Cl) 
2 
V--y-(2-p) ti 
ae Coe 
Oe M “2(atp-2y,) C722) 
Wie Gemmn Sino = | Lait 4p - VY(1+4p)7-1 


Equations (7.1) and (7.2) implies that x, and ome depends 
only on the ratio of A and uw, when M is kept constant. To 
check this implication, Simulation Program II was run with 


five sets of XA and yw where the ratio is kept equal to 2. 


as 
2) 
The steady state means and variances obtained from these 

runs turned out to be extremely close to x and ce as calcu-~ 


lated from equations (7.1) and (7.2) and these results are 


given in Table I. 


Sil 





Table I 


Steady State Mean and Variance of X(t) with 


o constant, N = 100, p = mn = 2.0. 


Arrival rate | Departure rate Simulation Diffusion 
of messages of messages "loss" model "loss' model 


A I= x Ce Xo oe 
(# per hour) (# per hour) z 
: : 13.369 : 13.389 ‘ 





(Each simulation runs for 2400 hours) 
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Next, keeping M=10, or N=20, Six more runs were made with 
the values of 9 ranging from 0.5 to 10. Again Xo and O° 
calculated from equations (7.1) and (7.2) were in close agree- 
ment with the simulation results. These results are summa- 
rized in Table II. 

Equations (7.1) and (7.2) also imply that Xo and oe are 
linear functions of N keeping p constant at 2. Six simula- 
tion runs were made for values of N ranging from 10 to 200. 
The steady state probabilities for the birth and death "loss" 
model were also calculated for the same values of pop and N. 
Again, the results indicated that all three models agree very 
well as far as steady state mean, variance and cumulative 
probabilities are concerned. The steady state means and 
variances for the three models are given in Tables III and 
IV respectively. In Table V we compared the cumulative 
probabilities of the three models for the case where N=20 
and p=2. These same values when plotted on normal probabil- 
ity graph paper (Figure 3), gave clear visual evidence of 


their close agreement. 
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Fable it 


Steady State Mean and Variance of X(t) with 


N constant, N = 40, A = 2.5. 


Departure rate Simulation Diffusion 
"loss" model . "loss" model 


‘ 
S S S S 
: 3.647 10.000 : 





(Each simulation runs for 2400 hours) 
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Table III 


Steady State Mean of X(t) with A and u constant, 


A = 2.5, p = 10.0 and p = i= = 2.0. 













Steady State Mean, Xo 
Number of 
mere Simulation Diffusion Birth and Death 
"loss" model | "loss" model "loss" model 








(Each simulation runs for 2400 hours) 
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Table IV 


Steady State Variance of X(t) with A and i 


COnseant, A = 2.50, 2 = 10.0 and 9» = oy = 2.0. 


WW) 











Steady State Variance o” 
Number of : 
ees Simulation | Diffusion Birth and Death 
"loss" model | ''loss" model "loss'' model 


1,559 
3.106 
.693 


elog 


1 Oo © 


. 746 


15.484 


(Each simulation runs for 2400 hours) 
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Table V 


Steady State Cumulative Probabilities, N=20, 


M=N/2=10, A = 2.5, p = 10.0, p = 4 = 2.0. 


Steady State Cumulative Probabilities P. 


Simulation Diffusion Birth and Death 
"loss" model | "loss" model "loss" model 





(Each simulation runs for 2400 hours) 
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VIII. APPROXIMATION OF "NO LOSS' MODEL BY "LOSS" MODEL 


Since the diffusion "loss" model gives relatively simple 


formulas for x_, oe 


and ay it appeared desirable to re- 
late the results of the "loss" model to the "no loss" model. 
It may be noted that in the "loss" model, only a fraction of 
the arriving messages are transmitted, whereas in the GPSS 
‘no loss" model all arriving messages are eventually trans- 
mitted. 

First, we need to determine the fraction of messages 
transmitted in order to obtain the effective rate of message 
arrivals at a "loss" model. In the "loss" model, the prob- 
ability that a message arrives at a busy machine iS approxi- 
mately Xo divided by M or Vins assuming steady state. Since 
a message can be transmitted only if both the transmitting 
and receiving machines are available, the probability of a 
message being transmitted immediately is approximately 
CL assuming independence between availability of ma- 


chines and that N is large. Hence ro is the effective 


fet 
arrival rate of messages at a machine in the "loss" model is 
given by 


d ort = CSN (8.1) 


This leads to the idea that a 'no loss" model with input 
parameters N, A and yw can be approximated by a "loss" model 
with input parameters N, A' and uw where A' is such that 


NE | = 2. Using equation (4.15) for ee with p' Sane 
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the equation AGN) aye = X} can be rewritten as 


MO (GE ee *~1 - 7G ae (8.2) 


After more rearrangement of the various terms (details are 
given in Appendix E), we obtain an equation expressing i' 


explicitly as a function of » and i: 


Ne | (8.3) 


The above relationship between "loss" and "no loss" mod- 
els was tested using the simulation models with five differ- 
ent values of u and A' as input parameters, N and A being 
kept constant at 20 and 2.5 respectively. The steady state 
means and variances are summarized in Table VI. For the 
case of N=20, w=10.0, A=2.5 and A'=10.0, the steady state 
cumulative probabilities are also calculated using the dif- 
fusion "loss" model, and the results from all three models 
are tabulated in Table VII and plotted in Figure 4. All the 
results indicate that the approximation of the 'no loss" 


model by the "loss" model is quite accurate. 
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Table VI 


Comparison of "no loss" model with its 


"equivalent" "loss" model, N=20, A=2.5. 





(Each simulation runs for 2400 hours ) 
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Zable Vif 


Steady State Cumulative Probabilities Calculated 


from Three Different Models, N=20, u=10.0. 


Steady State Cumulative Probabilities ae 
State 


Simulation Simulation 
"no loss" model |} "loss" model 


A=2.5 A'=10.0 










Diffusion 
"loss" model 
A'=10.0 







. 0002 
. 0034 
.0266 
.1146 
. 3492 


.6507 


oOo. OC @ Oo 36.) 2 Se 


.8874 
.9734 


7966 


- 2 oO 0 0 0D OD Oo Oo © 


. 9998 






eS Oe Ono BO Ore ee So 
ne ee Se 


— 


. 0000 


(Each simulation runs for 2400 hours) 
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IX. CONCLUSION 


It was found that the steady state results of the diffu- 
sion and the birth and death models were in very close agree- 
ment with the simulation model, through comparison of the 
steady state means, variances and cumulative probabilities. 
The usefulness of the diffusion model, as compared to the 
birth and near model, lies in the ease with which its steady 
state mean and variance and its discretized probabilities 
can be computed without resorting to a digital computer. 

The simple form of equations (7.1) and (7.2) exhibits 
the fact that the steady state mean and variance vary line- 
arly with N and depend only on pop, and not on the individual 
values of A and uw. These results cannot be inferred from 
the steady state formulas of the birth and death model or 
the simulation model. 

Lastly, we developed a simulation model for the much 
more complicated case when the messages are allowed to wait 
for transmission. We showed that a "no loss" model with in- 
put parameter N, A and uw can be quite accurately approxi- 
mated by a "loss" model with input parameter N, A', and uu 
where A‘' is calculated through equation (8.3) from the other 


parameters. 
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APPENDIX A 


Derivation of the Characteristic Function of S(t) 


Define 4,08, t) and Y,, 08, t) as the characteristic func- 


tions of X(t) and 54,0) respectively: 


$,,(0,t) 


YO, t) 


Since 


10X,,(t) 
Ele ] 2 Ce ®) 


il 


10S, ,(t) 
Ele ] (A.2) 


ee) hiya aa /M Sy(t) 


we obtain 


Vy (Ot) = 


oe” OY My (t) o,(O/vM, t) (A.3) 


Given that X(t) = jJ, equations (4.3), (4.4) and (4.5) give 


the respective probabilities that X,,(ttdt ) be equal to jtl, 


j=2. (Or “Ge 


Hence we can write down the following expression 


for 9,00, t+dt): 


$y (8, t+dt) 


{I 


Bret SMC Eat), 


Fy HOCSCE AY 7: ave ddt_ 
E[e {4M?-2M+2X,,(t) (1-4M)+4X,,(t) ay 


+ e* PMC) (42x () pAdt~(4M2—2M+2X, (t) (1-4M) 


Adt 


2 3 
FAR Ct) ony 





1+ oF PANO I“T) fox (ty padt] + o(dt)] 
(A.4) 
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Then after rearranging the terms and taking the limit as 


t * ©, we obtain 


2b 4(8,t) am oy CO. ttat)-6,,(8, t) 
dt dt+0 dt 


= pretOXu(t) 16; 


4M* -2M+2X, (Ce) le 4M)+4x?2 2 (t) gy 


1 


+ oF PAMOT) 90x (t)- (4M? -2M+2X, (t) (1-4M)+4X2(t)) 





a 


LO Xa t ) - 20 
2M—1 . 


Further rearranging of terms results in 


3,,(0,t) 
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= 2 LOAn(t),.10 ,, 4A _ 
at E[X,, (te (e 1)SM-1 
+ 2rX,(t o> AMC) {SESE ae (e pom) 
a 2,merOXm(t) (48 _ 4} 
However, Since 
oo,,(6,T) ; 
M ‘ = ‘ 1OXy(t) 
as ae E[ix,,(t)e J 
aad co, tf) 
M*~? 1O0Xy(t) 
= = M 
nye ap X,, (te ] 
the above can now be written as 
OO > C ) eedee Const) ,; 
M M oi 8 4x 16 
ee Sa = 
at 762 ca oui — A 
06,,(6,t) 
hae) aid 4M-1 -16 


CAD) 


CAG) 


CAz 





Rewriting equation (A.3) we obtain the following expressions 


for o, (8 /YM, t) and its first and second partial derivatives 


oibYMy(t), 


,(8/YM,t) = 60> t) (A.8) 


sce aa 


at + i0VMy'(t)¥,,(6,t)} (A.9) 


96,,(8//M, t) : ay.,(9,t) 
— — = Fo aged) ea Maa + ivfly(t)¥,,(0,t)} (A.10) 
07o,,(6//M, t) 2) a (Ce 1) aV,,(68,t) 
557 Peres BT Ft 2iviy(t) x6 
—My*(t)p,,(6, t)3 (A119) 


If we rewrite equation (A.7) with 6/YM in place of 8, we 


obtain 


at e 62 1- 


Ody (O/YM,t) — 8° dy (8/VM, t) 910 /Vi 4am 


2M-—1 


oy COME) an i0/VM_, 


~i0//M 
v6 Bi-i‘° : “aes 


+ 2/Mi )-2¢ 


i6//M_, 


+ 22M, (O/VM, t)(e ) (CAael») 


Using equations (A.8) through (A.11) to replace all the terms 


i6/My(t) 


involving > by wW, and factoring out e equation 


(A.12) becomes 


47 





aU (Oy ae 
———— + i6/My'(t)W,,(6,t) 
4M 


aoe 2 
a io /V/M a Vy 69 >t) UC 
~ 2M-1 


) ea + 2i/My(t) 





(l-e —My* (tp, (6 Sb) 


06 


~io/yM., ,°¥m 69» *) 


4M-1, i0//M 
re 


OMc 1S © —-1)+p (1-e 


+ 


2r/Mif 


+ivliy(t)yy(o.t)) + 2rmce e/a yy (0, t) (A.13) 


substituting eo iva and gee il with their power series 


= 2 
Se ee OG i 


5 - elb/ vm 
13/2 





x 


1 en ie/vM _ i6 
v7 43/2 


YM 


+2_ + O( 





and after more rearranging of terms in (A.13) we may write 


ae io Vily fy! (t)-2a (sey? (t)- (Sas +o )y(t)-1)} 





2 
4\M .-ie , 02.°%M 4M 





162 
rae Gua 7) nee e Sa {os-zy(t)(o+ =+) 
2M-1 JM 2M’ 96 2M-1 2 /M 
- (55) (0+ 107 )45¢-94 182yy0, 0M _ (aati y?(t) 
a Ji oi oe: 
= 2 
- (5 “p)¥CLIFLIAND- yy) (A.14) 


where the arguments of Wy,0 0, ) have been dropped for ease of 
exposition. The deterministic process y(t) must satisfy 


the following first order differential equation, as we let 


M+ © 


Wea crey=(t)>— (2ro)y(t) -1) = 0 (A.15) 
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in order that equation (A.14) simplify to 


Bb = ~92¥for(atp-2y(t)}- Sp {2r(y2(t)-(2-p)y¥(t)-1)) 
(A.16) 
This is the transformed version of the forward differential 
equation for a non-Stationary O.U. process. 

We have to solve the first order differential equation 
(A.15) in order to obtain y(t) as explicit function of t. 
For the sake of convenience, y'(t) and y(t) are written as 
dy 


at and y respectively. Then equation (A.15) is 


CY = 2a{y?-y(2+p)+1} (A.17) 


Now let 
Vee lee p 
and 


w2 = (1439)? - 1 


Then equation (A.17) can be written as 


chy ree 
Ge = 2a{ (y-v)?-w?} 


which in turn may be written as 


1 t, 1 
y-V-w y-vtw 








( ydy = 4rAwdt 


Integrating both sides of the equation we find 


In 2 = anwt + c, 
y-Vv+tw 


or 


— y+ ce 
VAT ant 4rAwt 


y-V-w 
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Assume that the system starts off empty, hence y=0O when t=0, 


then c2 1s evaluated as 


mee. W-V 
: Wtv 
Now let 
pee wee Pei Re ey = VE 
Wtv y-V-w 
Then 
- 1+k 
y V- ick W 


Rewriting y as a function of t again 


7 ree iY -4,r>wt 
= iene WIV ean 
WiC a 7 ET Weer ee (A.18) 
7” wty 
Further, 
=. lim y(t) _ vy (A.19) 
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APPENDIX B 


Derivation of Stochastic Differential Equation 


From equations (4.21), (4.22) and (4.23) 


dX,,(t ) = Bebe tat ae Z(t )V¥a(x,t dt (Bei 
{2M-2X,,(t)}{2M-2X,,(t)-1} 

By (x,t) = =) ss a ooo ~2X,,(t)pr (B22) 
{2M-2X,,(t)}{2M-2X,,(t )-1} 

Oy CX, t) = SE ae 1 ee ae +2X,,(t pr CB ac) 


where Z(t) is a purely random Gaussian process with zero 
mean and unit variance. Using equation (4.1) we can express 
X(t) in terms of S46 t) 

X(t) = VMS, (t) + My(t) 
therefore 

shee) = YMdS,,(t ) + My'(t)dt (B.4) 
Using equations (B.2), (B.3) and (B.4) to substitute for 
Mew anc (x. cand dX, Ct) in equation (B.1) we obtain 
YMdS,,(t) + My'(t)dt 


{2M-2/MS,,(t )-2My(t)}{2M-2/MS,,(t)-2My(t)-1} 


= i aaa hdt 


2pr{ YMS,,(t )+My(t) }dt} + Z(t )[2patvMS,,(t)+My(t) }dt 


{2M-2/MS,,(t)-2My(t)}{2M-2/MS,,(t )-2My(t)-1)} 


4 
2M-1 Adtj~ 


at 





Carrying out all the necessary multiplication, rearranging 
some of the terms and (Since we are eventually going to let 
M go to infinity) approximating 2M-1 by 2M, we obtain the 
following expression that is exact in the limit. 


YMdS,,(t) + My'(t)dt 


ons 


a iL 
= Adt{ 2M-4My (t)+2My* (t)+S,,(t ){4/M(y(t)-1) + 
+28,,7(t)-20{/MS,,(t )+My(t)}] + Z2(t)VrAdtl2M-4My(t )+2My?(t) 


a 
2 


= Il 7 
+S (4) {4/M (y(t )-1) + eh+28,° (t)+2p{ VMS, (t )+My(t ) 3] 


Dividing by /M throughout the whole equation and after more 


rearranging of terms, we have 


dS, (t) + /My'(t)dt 


2 


aviirat{y? (t)-(2+p y(t) +1} +28, (t)dt{2y(t)-2-p+ ay 


2X 


Ty Sy (t) + Z(t) 2rdtfy?(t)-(2-p y(t) +1} 


“te 


2\5,,(t )dt | 
M 1 2h 


ae 
2 


Now we let M go to infinity and assume that y'(t) = 2Afy?(t) 


(t) we have. 


-(2+p)y(t)t+1}. Then with S(t) = lim Sy 


M>© 
d5(t) >= 2d{2y(t)-2-p}S(t)dt+Z(t)[2rLy2(t)-(2-p) y(t +1} dt]? 


Again we have shown that the process S(t) is an Ornstein- 


Uhlenbeck process with the following infinitesimal mean and 
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variance, conditioned on S(t)=s. 


-B(t)s = -2A{2y(t)-2-p}s 


BCs, t) 


a(S,t) = a(t) = 2Aly?(t)-y(t)(2-0)+1} 
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APPENDIX C 


Derivation of Steady State Mean and Variance of X(t) 


From equation (4.6) we have 


2 
St = _a(t)es¥ — a(t)(S-)v (C.1) 
where 
B(t) = 2A{2+p-2y(t)} 


a(t) = 2r0f{y*(t)-y(t)(2-p)+1} 


Let the mean and variance of S(t) be denoted by u(t) and 
1(t), and since S(t) is normally distributed, we can write 
icmendibacterustie function W(0,t) in terms of W(t) and 


EC) 


_tidu(t)-20*1(t)} 


y(e8,t) = (C.2) 


Differentiating w(8,t) with respect to 9@ and t then dropping 


the arguments of w(6,t), we obtain 


BY = {iey'(t)-g02r' (t)}e! HOH(E) “20° 1(t)} (C.3) 
° z 
st = {in(t)-07(t)}et tou (t)-28°t(t)} (C.4) 


We substitute these equations for y, oe, oe into equation 


Bale 
(Cp ec ainceliangroutst he factor ol i6u(t)-20°rt(t)} 


we ob- 
ieee nl 


i6u'(t)-4071'(t) = -B(t)e{in(t)—-61(t)}-46°%a(t) (Ge5) 
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By equating coefficients of i8 and 6” on both sides, we ob- 


tain the following two first order differential equations: 


ict) = =Peepuct) . (C.6) 


Cle —28Ct 1 (to a(t) (C7) 


Equation (C.6) can be easily solved by separation of vari- 


ables [Ref. 7], yielding 


_j B(x)dx 
wis) =P Cee (C.8) 
It is reasonable to assume that the noise process S(t) starts 
off at time zero with zero mean and variance, hence y(0)=0 
and 1+(0)=0. 
Substituting the initial condition »(0)=0 into equation 


(C.8), we obtain 
u(t) = O noe At 0 (C.9) 


Equation (C.7) can be solved [Ref. 7] by the use of an inte- 
grating factor, giving us the following expression for T(t), 


with t(0)=0 


t Z 
=—2 VP bOx dx =v 2 f B(x)dx 
t(t) =e “ tf aczye f° dz } (C.10) 


Note that 8(t) and a(t) are functions of y(t) where y(t) is 


a complex function of t, given by equation (A.18): 


eas oe AAwt 


y(t) =v-wf{ fee (C.11) 


w-v _-—-4Awt 
e 


sy) 





Therefore the two integrals on the right hand side of equa- 
tion (C.10) are quite intractable, unless t is very large, 
in which case we are able to derive an approximation for 
mC t). 

For the typical values of A and p that we are using, the 
term Se in equation (C.11) vanishes very rapidly with t 
approaching infinity. For example, with A=2.5 messages/hour, 


a fawt . .-4(2.5)(73) (1) 


p=2eend t=1,0 hour: =<3 ns 0- 7) 


So for time t greater than one hour, y(t) can be very accu- 


rately approximated by ve where 
y, = 1+3p - V(1+3p)*-1 


and hence @(t) and a(t) can be approximated by the following 


expressions. 
B(t) = 2d (2tp-2y_,) = Be CC) 
eS 2 Bs = 
a(t) = 2X(¥, y (2 p)e®) a CCaLSy 


Referring to Figure 5, we can see that when z is very large, 
such that the area under the curve 8(x) up to time x=1.0 
hour, iS insignificant when compared to the area under the 
curve up to time x=z, then B 2 is a good approximation for 


Z 


J B(x)dx, 





(Anoy )X 


(x9 


7702S ox off 


O¢ 


G eansty 


(x)¢ 


0! 


‘SA X JO 40Td 


(ayy 


O7 





Z 
23 8 (x )dx 2B 2 
Similarly a(z)e ° can be approximated by ae 


Thus an approximate solution for T(t) is 


Oe -2B t 
Tu a= 5B. fl -e ] t > 10 hours 
eB 
2B t 
However, for such large value of t, e ~ 0, so we have 
Ce 
i) 2 ee t > 10 hours 
2B, = 


a — 
7 oe y (2 p)+1 
2(2tp-2y.) 


“~ 
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TSTOP 


APPENDIX D 


Notations used in Simulation Program I 


arrival time of next message at machine r 

departure time of next message at machine r, 

8.0 when machine r is not in use 

denotes the machine to which machine r is transmitting 
number of messages undergoing transmission at epoch qT; 
current clock time of simulation model 

exponential variate with mean = 1 

2uniVval rate Ol Messases AL each machine 

departure rate of a message undergoing transmission 
total number of machines 

discrete uniform (0,N) variate 


length of one simulation run 


o9 














a =minfa,} 


k = 
S =min{S, 8,70} 





CLOCK= a 


=a + 
a, =a, e/ir 








Generate 
n, n#k 





= CLOCK + e/u 


L 
L 
S 
S = CLOCK + e/u 


k 
n 
k 
n 
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Notations used in Simulation Program II 


a arrival rate of messages at each machine 

d departure rate of a message undergoing transmission 

N total number of machines 

M = N/2 

ds | rate of entering state i 

ie rate of leaving state i 

TST initial state at time 0 

N. number of messages undergoing transmission at epoch qT; 
e exponential variate with mean = 1 

| uniform (0,1) variate 


TSTOP length of one simulation run 
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Listing of GPSS Program 


This program simulates the arrivals and departures of 
messages at a system of 20 sendox machines. Functions EXP 
generates exponential variates with mean of one and UNI gen- 
erates discrete uniform variates between 1 and 20. The pro- 
gram consists of 20 closely similar segments. freer segment 
Simulates the activities at one machine. At the end of each 
Simulation run, this program gives the expected waiting 
times at each machine and steady state probabilities of 


X(t) - 
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APPENDIX E 


Calculation of A' 


In the previous sections, we had shown that a "no loss" 
model with input parameters N, A and u can be approximated 
by a "loss" model with input parameters N, A' and u, where 
wes chosen such that the effective arrival rate CN) 
for the "loss" model, iS equal to id. 

When X(t) is in steady state, the mean number of busy 
machines 1s Ny.» hence the probability of finding a busy 
machine is Ny ,/N or y.- In the "loss" model, a message is 
transmitted only if it finds both the transmitting and re- 
ceiving machines free. .The probability of this event is 
approximately Gey) assuming independence between avail- 


ability of machines and a large N. Therefore, the effective 


Bea ay 2 rate, A pep, for a "loss" model is Giey = where 
We 1s given by equation (4.15) as 
yi, = ithe! - /(1+8—)?-1 
fier fac 
with op W. 


Therefore, we require that 
FOGG Moyo le 3 oN) ae ee 


Poreweec: pandine the quadratic term in equation (E.1) and re- 


arranging terms, we obtain 





(1+4e') - shor = Ye'CHe) (E. 2) 
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Squaring both sides of equation (E.2) and recalling that 


o= sf we obtain 
1 
a a 
p p ; 


Making the substitution for p and o', we finally have 


Using this rate of arrivals in our "loss' model we are able 
to approximate the system size characteristics of the "no 


loss" model with arrival rate i. 
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